Summarise Substrate Data

Begin

Read in libraries and 2016-2019 data; 2015 UPC uses quadrats and is inconsistent with later methods.

#to load data make sure directory is set to knit from project directory
knitr::opts_chunk$set(message=FALSE, warning=FALSE)

# Libraries
library(tidyverse)
library(ggplot2)
library(viridis)
library(reshape2)
library(janitor)
library(magrittr)

#base.dir <- "/Users/ole.shelton/GitHub/OCNMS"
base.dir <- "/Users/jameal.samhouri/Dropbox/Projects/In progress/OCNMS/OCNMS"
data.dir <- paste0(base.dir,"/Data/CSV_2015_on")
data.summary.output.dir <- paste0(base.dir,"/Data/Summarized data")
setwd(data.dir)

# 2015 data. tricky because in 2015 we took % cover data in quadrats, and that included things that we/PISCO now call substrate but also that we/PISCO now call % cover. exclude for now
# dat_2015 <- clean_names(read.csv("2015_OCNMSDataComplete_standardized_122116.csv"), case = "all_caps") %>%
#   remove_empty(c("rows", "cols")) 
# glimpse(dat_2015)
# head(dat_2015)
# upc_2015 <- dat_2015 %>% filter(PISCO_DATATYPE=="UPC")
# head(upc_2015)


dat.2016.on.upc <- clean_names(read.csv("NWFSC_UPC_ALLYEARS_data_2019.csv"), case = "all_caps") %>%
  remove_empty(c("rows", "cols")) 
# fix data entry issues
dat.2016.on.upc %<>% mutate(
  CLASSCODE = toupper(CLASSCODE)
)
unique(dat.2016.on.upc$CLASSCODE)

##  [1] "BEDRK"     "BOULD"     "SAND"      "0-10CM"    "10CM-1M"   "BARSAN"   
##  [7] "LEAFYRED"  "CRUCOR"    "ZOSTERA"   "BROWNUN"   "PTERO"     "COB"      
## [13] "MACPYR_HF" "HYDROID"   "NEREO_HF"  "1M-2M"     "LACY"      "SHELL"    
## [19] "BRANCH"    "BARROC"    "ARTCOR"    "CUPCOR"    "BARNAC"    "BRYO"     
## [25] ">2M"       "SPONGE"    "SOLTUN"    "BIVALVE"   "LAMHOLD"   "PHYSPP"   
## [31] "SCUM"      "COLTUN"    "BUSHY"     "TUBEWORM"  "ENCRED"    "SCALLOP"  
## [37] "URCHIN"    "HYDROCOR"  "GREEN"     "CUCSPP"    "ANEM"

unique(dat.2016.on.upc$SITE)

## [1] Neah Bay           Destruction Island Cape Johnson       Tatoosh Island    
## [5] Cape Alava        
## 5 Levels: Cape Alava Cape Johnson Destruction Island ... Tatoosh Island

glimpse(dat.2016.on.upc)

## Rows: 6,570
## Columns: 17
## $ ENTRY_ORDER <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17…
## $ YEAR        <int> 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 201…
## $ MONTH       <int> 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, …
## $ DAY         <int> 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 1…
## $ SITE        <fct> Neah Bay, Neah Bay, Neah Bay, Neah Bay, Neah Bay, Neah Ba…
## $ OBSERVER    <fct> Kinsey Frick, Kinsey Frick, Kinsey Frick, Kinsey Frick, K…
## $ BUDDY       <fct> Steve Lonhart, Steve Lonhart, Steve Lonhart, Steve Lonhar…
## $ SIDE        <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ ZONE        <int> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, …
## $ TRANSECT    <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ HEADING     <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ DEPTH_FT    <int> 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 1…
## $ SEGMENT     <fct> 0-10m, 0-10m, 0-10m, 0-10m, 0-10m, 0-10m, 0-10m, 0-10m, 0…
## $ CATEGORY    <fct> SUBSTRATE, SUBSTRATE, SUBSTRATE, RELIEF, RELIEF, COVER, C…
## $ CLASSCODE   <chr> "BEDRK", "BOULD", "SAND", "0-10CM", "10CM-1M", "BARSAN", …
## $ COUNT       <int> 2, 2, 6, 9, 1, 5, 3, 1, 1, 6, 4, 4, 6, 3, 3, 2, 1, 1, 2, …
## $ NOTES       <fct> , , , , , , , , , , , , , , , , , , , , , , , , ,

substrate.codes <- clean_names(read.csv("substrate_codes.csv"), case = "all_caps") %>%
  remove_empty(c("rows", "cols")) # js checked on 072720 and codes are same in 2019 data
  
substrate <- dat.2016.on.upc %>% filter(CATEGORY=="SUBSTRATE")
#relief <- dat.2016.on.upc %>% filter(CATEGORY=="RELIEF")
#cover <- dat.2016.on.upc %>% filter(CATEGORY=="COVER")

Analyze substrate data

Summarise raw data into data frames that: 1) sum segments for each transect. 2) calculate summary stats for site-year-zone. 3) calculate summary stats for site-year. 4) calculate summary stats for site-zone. 5) calculate summary stats for site.

# Pad with zero observations. i.e., complete the data so that each segment has all 4 substrate categories
  cat <- data.frame(merge(unique(substrate$SITE),unique(substrate$CLASSCODE))); colnames(cat) <- c("SITE","CLASSCODE")
  dat.sub <- substrate %>% 
    group_by(YEAR,SITE,SEGMENT,TRANSECT,SIDE,ZONE,OBSERVER) %>% 
    summarise(a=length(SEGMENT)) %>% 
    full_join(.,cat,by="SITE") %>% 
    dplyr::select(-a) %>% 
    left_join(.,substrate)
  dat.sub$COUNT[is.na(dat.sub$COUNT)==TRUE] <- 0
  
  # order sites from south to north
  dat.sub$SITE <- factor(dat.sub$SITE,
               levels=c("Destruction Island","Cape Johnson","Cape Alava","Tatoosh Island","Neah Bay"))
  
  # make year a factor
  dat.sub$YEAR <- as.factor(dat.sub$YEAR)
  
  # not sure if we have to account for SIDE 
  # (eg, multiple observers on a segment)? looks like no.
  tmp<-substrate %>% 
    group_by(YEAR,SITE,SEGMENT,TRANSECT,SIDE,ZONE,OBSERVER) %>% 
    summarise(a=length(SEGMENT), b=length(OBSERVER))
  which(tmp$a != tmp$b) # integer(0)

## integer(0)

  # with(dat.sub, table(SITE, ZONE, TRANSECT, SIDE))
  # unique(dat.sub$SIDE)
  # length(which(dat.sub$SIDE != 1))/dim(dat.sub)[1]
  
  # but it does look like there are some issues with duplicate entries. so fix those
  dat.sub <- dat.sub %>% 
    group_by(YEAR,SITE,ZONE,SIDE,TRANSECT) %>%
    distinct(CLASSCODE, SEGMENT, COUNT, .keep_all = TRUE) %>% # remove duplicate entries based on CLASSCODE, SEGMENT, COUNT
    mutate(
      TOTAL_COUNT = sum(COUNT),
      PROPORTION = COUNT/TOTAL_COUNT
    )
  unique(dat.sub$TOTAL_COUNT) # 30 29

## [1] 30 29

  View(dat.sub[which(dat.sub$TOTAL_COUNT !=30),])
  glimpse(dat.sub)

## Rows: 3,324
## Columns: 19
## Groups: YEAR, SITE, ZONE, SIDE, TRANSECT [277]
## $ YEAR        <fct> 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 201…
## $ SITE        <fct> Cape Alava, Cape Alava, Cape Alava, Cape Alava, Cape Alav…
## $ SEGMENT     <fct> 0-10m, 0-10m, 0-10m, 0-10m, 0-10m, 0-10m, 0-10m, 0-10m, 0…
## $ TRANSECT    <int> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, …
## $ SIDE        <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ ZONE        <int> 5, 5, 5, 5, 10, 10, 10, 10, 5, 5, 5, 5, 10, 10, 10, 10, 5…
## $ OBSERVER    <fct> Kelly Andrews, Kelly Andrews, Kelly Andrews, Kelly Andrew…
## $ CLASSCODE   <chr> "BEDRK", "BOULD", "SAND", "COB", "BEDRK", "BOULD", "SAND"…
## $ ENTRY_ORDER <int> NA, 1037, NA, 1038, 883, 884, 886, 885, 1058, 1059, NA, 1…
## $ MONTH       <int> NA, 8, NA, 8, 8, 8, 8, 8, 8, 8, NA, 8, 8, 8, 8, 8, NA, 8,…
## $ DAY         <int> NA, 10, NA, 10, 10, 10, 10, 10, 10, 10, NA, 10, 10, 10, 1…
## $ BUDDY       <fct> NA, Steve Lonhart, NA, Steve Lonhart, Greg Williams, Greg…
## $ HEADING     <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ DEPTH_FT    <int> NA, 22, NA, 22, 31, 31, 31, 31, 16, 16, NA, 16, 26, 26, 2…
## $ CATEGORY    <fct> NA, SUBSTRATE, NA, SUBSTRATE, SUBSTRATE, SUBSTRATE, SUBST…
## $ COUNT       <dbl> 0, 7, 0, 3, 2, 5, 1, 2, 1, 7, 0, 2, 7, 1, 1, 1, 0, 8, 0, …
## $ NOTES       <fct> NA, , NA, , , , , , , , NA, , , , , , NA, , NA, , , , , ,…
## $ TOTAL_COUNT <dbl> 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 3…
## $ PROPORTION  <dbl> 0.00000000, 0.23333333, 0.00000000, 0.10000000, 0.0666666…

  # create summary df. 1) sum segments for each transect. 2) calculate summary stats for site-year-zone. 3) calculate summary stats for site-year. 4) calculate summary stats for site-zone. 5) calculate summary stats for site
  #proportion data so use binomial distribution for summary stats. https://sites.warnercnr.colostate.edu/gwhite/wp-content/uploads/sites/73/2017/04/BinomialDistribution.pdf
  
  # 1) sum segments for each transect to get total proportion of each substrate CLASSCODE
  sub.summary.year.transect <- dat.sub %>% 
    group_by(YEAR,SITE,ZONE,SIDE,TRANSECT,CLASSCODE) %>%
    dplyr::summarise(
      PROPORTION=sum(PROPORTION),
      .groups = "keep"
    )
  sub.summary.year.transect$YEAR <- as.factor(sub.summary.year.transect$YEAR)
  # order sites from south to north
  sub.summary.year.transect$SITE <- factor(sub.summary.year.transect$SITE,
                                       levels=c("Destruction Island","Cape Johnson","Cape Alava","Tatoosh Island","Neah Bay"))
  glimpse(sub.summary.year.transect)

## Rows: 1,108
## Columns: 7
## Groups: YEAR, SITE, ZONE, SIDE, TRANSECT, CLASSCODE [1,108]
## $ YEAR       <fct> 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016…
## $ SITE       <fct> Destruction Island, Destruction Island, Destruction Island…
## $ ZONE       <int> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 10, 10, 10…
## $ SIDE       <fct> 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1…
## $ TRANSECT   <int> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 1, 1, 1, 1…
## $ CLASSCODE  <chr> "BEDRK", "BOULD", "COB", "SAND", "BEDRK", "BOULD", "COB", …
## $ PROPORTION <dbl> 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0…

  # 2) calculate summary stats for site-year-zone
  sub.summary.year.zone <- sub.summary.year.transect %>% 
    group_by(YEAR,SITE,ZONE,CLASSCODE) %>%
    dplyr::summarise(
      MEAN=mean(PROPORTION),
      SD=sqrt(MEAN * (1 - MEAN)),
      N=length(PROPORTION),
      SE=SD/sqrt(N),
      .groups = "keep"
      )
  sub.summary.year.zone$YEAR <- as.factor(sub.summary.year.zone$YEAR)
  # order sites from south to north
  sub.summary.year.zone$SITE <- factor(sub.summary.year.zone$SITE,
                         levels=c("Destruction Island","Cape Johnson","Cape Alava","Tatoosh Island","Neah Bay"))
  glimpse(sub.summary.year.zone)

## Rows: 156
## Columns: 8
## Groups: YEAR, SITE, ZONE, CLASSCODE [156]
## $ YEAR      <fct> 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016,…
## $ SITE      <fct> Destruction Island, Destruction Island, Destruction Island,…
## $ ZONE      <int> 5, 5, 5, 5, 10, 10, 10, 10, 5, 5, 5, 5, 10, 10, 10, 10, 5, …
## $ CLASSCODE <chr> "BEDRK", "BOULD", "COB", "SAND", "BEDRK", "BOULD", "COB", "…
## $ MEAN      <dbl> 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.91666667,…
## $ SD        <dbl> 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.2763854, 0.19…
## $ N         <int> 4, 4, 4, 4, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,…
## $ SE        <dbl> 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.11283387,…

  # 3) calculate summary stats for site-year
  sub.summary.year <- sub.summary.year.transect %>% 
    group_by(YEAR,SITE,CLASSCODE) %>%
    dplyr::summarise(
      MEAN=mean(PROPORTION),
      SD=sqrt(MEAN * (1 - MEAN)),
      N=length(PROPORTION),
      SE=SD/sqrt(N),
      .groups = "keep"
    )
  sub.summary.year$YEAR <- as.factor(sub.summary.year$YEAR)
  sub.summary.year$SITE <- factor(sub.summary.year$SITE,
                                       levels=c("Destruction Island","Cape Johnson","Cape Alava","Tatoosh Island","Neah Bay"))
  glimpse(sub.summary.year)

## Rows: 80
## Columns: 7
## Groups: YEAR, SITE, CLASSCODE [80]
## $ YEAR      <fct> 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016,…
## $ SITE      <fct> Destruction Island, Destruction Island, Destruction Island,…
## $ CLASSCODE <chr> "BEDRK", "BOULD", "COB", "SAND", "BEDRK", "BOULD", "COB", "…
## $ MEAN      <dbl> 0.950000000, 0.023333333, 0.013333333, 0.013333333, 0.22777…
## $ SD        <dbl> 0.2179449, 0.1509599, 0.1146977, 0.1146977, 0.4193985, 0.49…
## $ N         <int> 10, 10, 10, 10, 12, 12, 12, 12, 12, 12, 12, 12, 8, 8, 8, 8,…
## $ SE        <dbl> 0.06892024, 0.04773771, 0.03627059, 0.03627059, 0.12106990,…

  # 4) calculate summary stats for site-zone
  sub.summary.zone <- sub.summary.year.transect %>% 
    group_by(SITE,ZONE,CLASSCODE) %>%
    dplyr::summarise(
      MEAN=mean(PROPORTION),
      SD=sqrt(MEAN * (1 - MEAN)),
      N=length(PROPORTION),
      SE=SD/sqrt(N),
      .groups = "keep"
    )
  sub.summary.zone$SITE <- factor(sub.summary.zone$SITE,
                                  levels=c("Destruction Island","Cape Johnson","Cape Alava","Tatoosh Island","Neah Bay"))
  glimpse(sub.summary.zone)

## Rows: 40
## Columns: 7
## Groups: SITE, ZONE, CLASSCODE [40]
## $ SITE      <fct> Destruction Island, Destruction Island, Destruction Island,…
## $ ZONE      <int> 5, 5, 5, 5, 10, 10, 10, 10, 5, 5, 5, 5, 10, 10, 10, 10, 5, …
## $ CLASSCODE <chr> "BEDRK", "BOULD", "COB", "SAND", "BEDRK", "BOULD", "COB", "…
## $ MEAN      <dbl> 0.978571429, 0.016666667, 0.003571429, 0.001190476, 0.92380…
## $ SD        <dbl> 0.14480811, 0.12801910, 0.05965462, 0.03448273, 0.26530263,…
## $ N         <int> 28, 28, 28, 28, 21, 21, 21, 21, 30, 30, 30, 30, 28, 28, 28,…
## $ SE        <dbl> 0.027366160, 0.024193335, 0.011273663, 0.006516624, 0.05789…

  # 5) calculate summary stats for site
  sub.summary.site <- sub.summary.year.transect %>% 
    group_by(SITE,CLASSCODE) %>%
    dplyr::summarise(
      MEAN=mean(PROPORTION),
      SD=sqrt(MEAN * (1 - MEAN)),
      N=length(PROPORTION),
      SE=SD/sqrt(N),
      .groups = "keep"
    )
  sub.summary.site$SITE <- factor(sub.summary.site$SITE,
                                  levels=c("Destruction Island","Cape Johnson","Cape Alava","Tatoosh Island","Neah Bay"))
  glimpse(sub.summary.site)

## Rows: 20
## Columns: 6
## Groups: SITE, CLASSCODE [20]
## $ SITE      <fct> Destruction Island, Destruction Island, Destruction Island,…
## $ CLASSCODE <chr> "BEDRK", "BOULD", "COB", "SAND", "BEDRK", "BOULD", "COB", "…
## $ MEAN      <dbl> 0.955102041, 0.035374150, 0.006122449, 0.003401361, 0.09712…
## $ SD        <dbl> 0.20708001, 0.18472363, 0.07800618, 0.05822191, 0.29612986,…
## $ N         <int> 49, 49, 49, 49, 58, 58, 58, 58, 64, 64, 64, 64, 50, 50, 50,…
## $ SE        <dbl> 0.029582859, 0.026389090, 0.011143740, 0.008317416, 0.03888…

  #### END BASIC SUBSTRATE SUMMARIES.

Write out the data frames we just made.

### WRITE OUT SUBSTRATE DATA FRAMES
  
  
  write_csv(sub.summary.year.transect, paste0(data.summary.output.dir, "/SUBSTRATE summary by YEAR-SITE-DEPTH-SIDE-TRANSECT 2016-2019.csv"))
  write_rds(sub.summary.year.transect, paste0(data.summary.output.dir, "/SUBSTRATE summary by YEAR-SITE-DEPTH-SIDE-TRANSECT 2016-2019.rds"))
  write_csv(sub.summary.year.zone, paste0(data.summary.output.dir, "/SUBSTRATE summary by YEAR-SITE-DEPTH 2016-2019.csv"))
  write_rds(sub.summary.year.zone, paste0(data.summary.output.dir, "/SUBSTRATE summary by YEAR-SITE-DEPTH 2016-2019.rds"))
  write_csv(sub.summary.year, paste0(data.summary.output.dir, "/SUBSTRATE summary by YEAR-SITE 2016-2019.csv"))
  write_rds(sub.summary.year, paste0(data.summary.output.dir, "/SUBSTRATE summary by YEAR-SITE 2016-2019.rds"))
  write_csv(sub.summary.zone, paste0(data.summary.output.dir, "/SUBSTRATE summary by SITE-DEPTH 2016-2019.csv"))
  write_rds(sub.summary.zone, paste0(data.summary.output.dir, "/SUBSTRATE summary by SITE-DEPTH 2016-2019.rds"))
  write_csv(sub.summary.site, paste0(data.summary.output.dir, "/SUBSTRATE summary by SITE 2016-2019.csv"))
  write_rds(sub.summary.site, paste0(data.summary.output.dir, "/SUBSTRATE summary by SITE 2016-2019.rds"))
  
  #### END WRITING OUT SUBSTRATE DATA FRAMES

MAKE SUBSTRATE PLOTS

1) plots for site-year-zone.

  # a) substrate type on x axis, site as facets
  sub.plot.year.zone_sub <- ggplot() +
    geom_jitter(data=sub.summary.year.transect,aes(y=PROPORTION,x=CLASSCODE,colour=YEAR),alpha=0.5, position=position_dodge(width=0.5))+
    geom_point(data=sub.summary.year.zone,aes(y=MEAN,x=CLASSCODE,fill=YEAR),color="black",size=3,shape=21, position=position_dodge(width=0.5))+
    geom_errorbar(data=sub.summary.year.zone,
                  aes(ymin=MEAN-SE-1e-10,ymax=MEAN+SE,x=CLASSCODE, colour=YEAR),width=0.1, position=position_dodge(width=0.5))  +
    scale_colour_viridis_d()+
    scale_fill_viridis_d()+
    facet_grid(SITE~ZONE) +
    theme_bw()
 plot(sub.plot.year.zone_sub)

  # b) site on x axis, substrate type as facets
  sub.plot.year.zone_site <- ggplot() +
    geom_jitter(data=sub.summary.year.transect,aes(y=PROPORTION,x=SITE, colour=YEAR),alpha=0.5, position=position_dodge(width=0.5) )+
    geom_point(data=sub.summary.year.zone,aes(y=MEAN,x=SITE,fill=YEAR),color="black",size=3,shape=21, position=position_dodge(width=0.5))+
    geom_errorbar(data=sub.summary.year.zone,
                  aes(ymin=MEAN-SE-1e-10,ymax=MEAN+SE,x=SITE, colour=YEAR),width=0.1, position=position_dodge(width=0.5))  +
    scale_colour_viridis_d()+
    scale_fill_viridis_d()+
    facet_grid(CLASSCODE~ZONE) +
    theme_bw() +
    theme(
      axis.text.x = element_text(angle = 45, hjust=1)
    )
  plot(sub.plot.year.zone_site)

  # c) compare zones, years as facets, lines connecting 5m v 10m
  
  sub.plot.year.zone_compare1 <- ggplot() +
    geom_jitter(data=sub.summary.year.transect,aes(y=PROPORTION,x=ZONE, colour=SITE),alpha=0.5, position=position_dodge(width=0.5))+
    geom_point(data=sub.summary.year.zone,aes(y=MEAN,x=ZONE,colour=SITE),size=2, position=position_dodge(width=0.5))+
    geom_line(data=sub.summary.year.zone,aes(y=MEAN,x=ZONE,colour=SITE),size=1)+
    geom_errorbar(data=sub.summary.year.zone,
                  aes(ymin=MEAN-SE-1e-10,ymax=MEAN+SE,x=ZONE,colour=SITE),width=0.1, position=position_dodge(width=0.5))  +
    scale_colour_viridis_d()+
    facet_grid(YEAR~CLASSCODE) +
    theme_bw()
  plot(sub.plot.year.zone_compare1)

  # d) compare zones, CLASSCODE as facets, no lines connecting 5m v 10m
  sub.plot.year.zone_compare2 <- ggplot() +
    geom_jitter(data=sub.summary.year.transect,aes(y=PROPORTION,x=ZONE,colour=YEAR),alpha=0.3 , position=position_dodge(width=0.5))+
    geom_point(data=sub.summary.year.zone,aes(y=MEAN,x=ZONE,fill=YEAR),size=2,shape=21, position=position_dodge(width=0.5))+
    geom_errorbar(data=sub.summary.year.zone,
                  aes(ymin=MEAN-SE-1e-10,ymax=MEAN+SE,x=ZONE,colour=YEAR),width=0.1, position=position_dodge(width=0.5))  +
    scale_colour_viridis_d() +
    scale_fill_viridis_d() +
    facet_grid(CLASSCODE~SITE) +
    theme_bw()
  plot(sub.plot.year.zone_compare2)

2) plots for site-year.

# a) substrate type on x axis, site as facets
  sub.plot.year_sub <- ggplot() +
    geom_point(data=sub.summary.year,aes(y=MEAN,x=CLASSCODE,fill=YEAR),colour="black",size=3,shape=21, position=position_dodge(width=0.5))+
    geom_errorbar(data=sub.summary.year,
                  aes(ymin=MEAN-SE-1e-10,ymax=MEAN+SE,x=CLASSCODE, colour=YEAR),width=0.1, position=position_dodge(width=0.5))  +
    scale_colour_viridis_d()+
    scale_fill_viridis_d()+
    facet_grid(SITE~.) +
    theme_bw()
  plot(sub.plot.year_sub)

  # b) site on x axis, substrate type as facets
  sub.plot.year_site <- ggplot() +
    geom_point(data=sub.summary.year,aes(y=MEAN,x=SITE,fill=YEAR),color="black",size=3,shape=21, position=position_dodge(width=0.5))+
    geom_errorbar(data=sub.summary.year,
                  aes(ymin=MEAN-SE-1e-10,ymax=MEAN+SE,x=SITE, colour=YEAR),width=0.1, position=position_dodge(width=0.5))  +
    scale_colour_viridis_d()+
    scale_fill_viridis_d()+
    facet_grid(CLASSCODE~.) +
    theme_bw() +
    theme(
      axis.text.x = element_text(angle = 45, hjust=1)
    )
  plot(sub.plot.year_site)

  # c) time series, CLASSCODE as facets
  
  sub.plot.year_ts1 <- ggplot() +
    geom_point(data=sub.summary.year,aes(y=MEAN,x=YEAR,colour=SITE, group=SITE),size=2, position=position_dodge(width=0.5))+
    geom_line(data=sub.summary.year,aes(y=MEAN,x=YEAR,colour=SITE, group=SITE),size=1, position=position_dodge(width=0.5))+
    geom_errorbar(data=sub.summary.year,
                  aes(ymin=MEAN-SE-1e-10,ymax=MEAN+SE,x=YEAR,colour=SITE, group=SITE),width=0.1, position=position_dodge(width=0.5))  +
    scale_colour_viridis_d()+
    facet_grid(.~CLASSCODE) +
    theme_bw()
  plot(sub.plot.year_ts1)

  # d) time series, SITE as facets
  
  sub.plot.year_ts2 <- ggplot() +
    geom_point(data=sub.summary.year,aes(y=MEAN,x=YEAR,colour=CLASSCODE, group=CLASSCODE),size=2, position=position_dodge(width=0.5))+
    geom_line(data=sub.summary.year,aes(y=MEAN,x=YEAR,colour=CLASSCODE, group=CLASSCODE),size=1, position=position_dodge(width=0.5))+
    geom_errorbar(data=sub.summary.year,
                  aes(ymin=MEAN-SE-1e-10,ymax=MEAN+SE,x=YEAR,colour=CLASSCODE, group=CLASSCODE),width=0.1, position=position_dodge(width=0.5))  +
    scale_colour_viridis_d()+
    facet_grid(.~SITE) +
    theme_bw()
  plot(sub.plot.year_ts2)

3) plots for site-zone.

  # a) substrate type on x axis, site as facets
  sub.plot.zone_sub <- ggplot() +
    geom_point(data=sub.summary.zone,aes(y=MEAN,x=CLASSCODE,fill=factor(ZONE)),colour="black",size=3,shape=21, position=position_dodge(width=0.5))+
    geom_errorbar(data=sub.summary.zone,
                  aes(ymin=MEAN-SE-1e-10,ymax=MEAN+SE,x=CLASSCODE, colour=factor(ZONE)),width=0.1, position=position_dodge(width=0.5))  +
    labs(fill='ZONE', colour='ZONE') +
    scale_colour_viridis_d()+
    scale_fill_viridis_d()+
    facet_grid(SITE~.) +
    theme_bw()
  plot(sub.plot.zone_sub)

  # b) site on x axis, substrate type as facets
  sub.plot.zone_site <- ggplot() +
    geom_point(data=sub.summary.zone,aes(y=MEAN,x=SITE,fill=factor(ZONE)),color="black",size=3,shape=21, position=position_dodge(width=0.5))+
    geom_errorbar(data=sub.summary.zone,
                  aes(ymin=MEAN-SE-1e-10,ymax=MEAN+SE,x=SITE, colour=factor(ZONE)),width=0.1, position=position_dodge(width=0.5))  +
    labs(fill='ZONE', colour='ZONE') +
    scale_colour_viridis_d()+
    scale_fill_viridis_d()+
    facet_grid(CLASSCODE~.) +
    theme_bw() +
    theme(
      axis.text.x = element_text(angle = 45, hjust=1)
    )
  plot(sub.plot.zone_site)

  # c) compare zones, years as facets, lines connecting 5m v 10m
  
  sub.plot.zone_compare1 <- ggplot() +
    geom_point(data=sub.summary.zone,aes(y=MEAN,x=factor(ZONE),colour=SITE, group=SITE),size=2, position=position_dodge(width=0.5))+
    geom_line(data=sub.summary.zone,aes(y=MEAN,x=factor(ZONE),colour=SITE, group=SITE),size=1, position=position_dodge(width=0.5))+
    geom_errorbar(data=sub.summary.zone,
                  aes(ymin=MEAN-SE-1e-10,ymax=MEAN+SE,x=factor(ZONE),colour=SITE, group=SITE),width=0.1, position=position_dodge(width=0.5))  +
    xlab("ZONE") +
    scale_colour_viridis_d()+
    facet_grid(.~CLASSCODE) +
    theme_bw()
  plot(sub.plot.zone_compare1)

  # d) compare zones, CLASSCODE as facets, no lines connecting 5m v 10m
  
  sub.plot.zone_compare2 <- ggplot() +
    geom_point(data=sub.summary.zone,aes(y=MEAN,x=factor(ZONE),colour=CLASSCODE, group=CLASSCODE),size=2, position=position_dodge(width=0.5))+
    geom_line(data=sub.summary.zone,aes(y=MEAN,x=factor(ZONE),colour=CLASSCODE, group=CLASSCODE),size=1, position=position_dodge(width=0.5))+
    geom_errorbar(data=sub.summary.zone,
                  aes(ymin=MEAN-SE-1e-10,ymax=MEAN+SE,x=factor(ZONE),colour=CLASSCODE, group=CLASSCODE),width=0.1, position=position_dodge(width=0.5))  +
    xlab("ZONE") +
    scale_colour_viridis_d()+
    facet_grid(.~SITE) +
    theme_bw()
  plot(sub.plot.zone_compare2)

4) plots for site

  # a) substrate type on x axis, site as facets
  sub.plot.site_sub <- ggplot() +
    geom_point(data=sub.summary.site,aes(y=MEAN,x=CLASSCODE, fill=CLASSCODE),colour="black",size=3,shape=21, position=position_dodge(width=0.5))+
    geom_errorbar(data=sub.summary.site,
                  aes(ymin=MEAN-SE-1e-10,ymax=MEAN+SE,x=CLASSCODE, colour=CLASSCODE),width=0.1, position=position_dodge(width=0.5))  +
    scale_colour_viridis_d()+
    scale_fill_viridis_d()+
    facet_grid(SITE~.) +
    theme_bw()
  plot(sub.plot.site_sub)

  # b) site on x axis, substrate type as facets
  sub.plot.site_site <- ggplot() +
    geom_point(data=sub.summary.site,aes(y=MEAN,x=SITE,fill=SITE),color="black",size=3,shape=21, position=position_dodge(width=0.5))+
    geom_errorbar(data=sub.summary.site,
                  aes(ymin=MEAN-SE-1e-10,ymax=MEAN+SE,x=SITE, colour=SITE),width=0.1, position=position_dodge(width=0.5))  +
    scale_colour_viridis_d()+
    scale_fill_viridis_d()+
    facet_grid(CLASSCODE~.) +
    theme_bw() +
    theme(
      axis.text.x = element_text(angle = 45, hjust=1)
    )
  plot(sub.plot.site_site)

PRINT OUT ALL OF THE PLOTS IN A SINGLE PDF.

  pdf(file=paste0(base.dir,"/Plots/Substrate v2.pdf"),onefile=T)
  
  print(sub.plot.site_sub)
  print(sub.plot.site_site)
  
  print(sub.plot.zone_sub)
  print(sub.plot.zone_site)
  print(sub.plot.zone_compare1)
  print(sub.plot.zone_compare2)
  
  print(sub.plot.year_sub)
  print(sub.plot.year_site)
  print(sub.plot.year_ts1)
  print(sub.plot.year_ts2)
  
  print(sub.plot.year.zone_sub)  
  print(sub.plot.year.zone_site)
  print(sub.plot.year.zone_compare1)
  print(sub.plot.year.zone_compare2)
  
  dev.off()

## quartz_off_screen 
##                 2